classdef simulator < handle
    properties
        N;                 % 总模拟光子数
        Ntr;               % 模拟的划分粒度
        light_source;      % 光源
        atmos;             % 大气
        K_sr;              % 瑞利散射系数，由于每个光子的波长不同，每个光子对应的散射系数都不同。
        K_sm;              % 米氏散射系数，由于每个光子的波长不同，每个光子对应的散射系数都不同。
        K_s;               % 散射系数 
        K_a;               % 吸收系数
        K_e;               % 衰减系数
        g;                 % 米氏散射非对称因子
        S1_vec;            % 米氏散射S1_vec
        S2_vec;            % 米氏散射S2_vec
        w_threshold;       % 光子湮灭能量阈值
        light_source_info; % 三维矩阵，前二维代表平面的x，y，第三维代表每个点的[x,y,光强，波长，方向余弦，电场相位]
        context;           % 三维矩阵，第一维代表第i个光子，第二维代表第j次散射，第三维是光子当前位置和状态[x,y,z,w,电场相位]  
    end

    methods(Static)
        % 得到方位角
        function angle = getAzimuth()
            angle = 2*pi*rand;
        end

        % 得到瑞利散射偏转角余弦值
        function cos_theta = getRandomRayleighDeflectionAngleCos()
            n = rand;
            cos_theta = nthroot((4*n-2)-sqrt((2-4*n)^2+1),3) + nthroot((4*n-2)+sqrt((2-4*n)^2+1),3);
        end
    end

    methods
        % 构造函数
        function obj = simulator(light_source, atmos)
            obj.N = 1000000;                                         % 500万个光子
            obj.Ntr = 1024;                                          % 模拟的划分粒度
            obj.light_source = light_source;                         % 光源
            obj.atmos = atmos;                                       % 大气
            obj.updateK(light_source.getLambda());                   % 散射和吸收，衰减系数初始化
            obj.g = 0.78;                                            % 米氏散射非对称因子（表征散射角余弦平均值，晴朗天气下的典型值0.78）
            obj.updateS(light_source.getLambda());                   % S1,S2数组初始化
            obj.w_threshold = 0.10;                                  % 光子湮灭能量阈值
            obj.light_source_info = zeros(18, obj.Ntr, obj.Ntr);     % 初始化光源信息三维矩阵 
            obj.context = zeros(4, 8, obj.N);                        % 初始化散射信息记录矩阵，最大散射次数为8次                     
        end

        % 更新对某个光子的散射系数
        function obj = updateK(obj, lambda)
            obj.K_sr = obj.atmos.get_K_sr(lambda);
            obj.K_sm = obj.atmos.get_K_sm(lambda);
            obj.K_s = obj.K_sr + obj.K_sm;
            obj.K_a = obj.atmos.get_K_a();
            obj.K_e = obj.K_a + obj.K_s;
        end

        % 更新对某个光子的S1,S2
        function obj = updateS(obj, lambda)
            radius = 0.05E-6;          % 50nm半径的微粒
            nmedium = 1.33;            % 介质折射率
            refrel = 1.5/nmedium;      % 有效折射率
            x = 2*pi*radius/lambda;    % size parameter
            nang = 90;
            [obj.S1_vec, obj.S2_vec, ~, ~, ~, ~] = mie(x, refrel, nang);
        end

        % 获取吸收自由程
        function s = getRandomDistance(obj)
            s = -log(rand)/obj.K_a;
        end

        % 获取散射自由程
        function s = getRandomStep(obj)
            s = -log(rand)/obj.K_s;
        end

        % 得到本次散射类型
        function scatter_type = getScatterType(obj)
            if rand <= obj.K_sr/(obj.K_sr + obj.K_sm)
                scatter_type = 0; % 瑞利散射
            else
                scatter_type = 1; % 米氏散射
            end
        end

        % 得到米氏散射偏转角余弦值
        function cos_theta = getRandomMieDeflectionAngleCos(obj)
            n = rand;
            cos_theta = (1+obj.g^2-((1-obj.g^2)/(1-obj.g+2*obj.g*n))^2)/(2*obj.g);
        end

        % 单个光子模拟
        function singleTimeSimulator(obj, index, lambda, x, y, z, u_x, u_y, u_z, m, n, s, E1, E2)
            % obj.updateK(lambda);
            % obj.updateS(lambda);
            k = 2*pi/lambda;
            w = 1.0;
            cnt = 1;
            % 获取吸收自由程，即最大传播距离
            dis = 0;
            distance = obj.getRandomDistance();
            single_context = zeros(4, 8);
            while w >= obj.w_threshold && dis <= distance && cnt <= 7
                single_context(:, cnt) = [x, y, z, E1];
                % 产生步长
                ss = obj.getRandomStep();
                % 更新坐标
                x = x + u_x*ss;
                y = y + u_y*ss;
                z = z + u_z*ss;
                % E1和E2跟随传播同步旋转
                E1 = E1*exp(1i*k*ss);
                E2 = E2*exp(1i*k*ss);
                dis = dis + ss;
                % 判断散射类型
                scatter_type = obj.getScatterType();
                % 产生新的方位角
                fi = obj.getAzimuth();
                % 产生新的偏转角余弦和正弦
                if scatter_type == 0
                    cos_theta = obj.getRandomRayleighDeflectionAngleCos();
                    sin_theta = (1-cos_theta^2)^0.5;
                else
                    cos_theta = obj.getRandomMieDeflectionAngleCos();
                    sin_theta = (1-cos_theta^2)^0.5;
                end
                % 更新方向余弦
                u_threshold = 0.9999999;
                if u_z <= u_threshold
                    u_x_new = sin_theta*(u_x*u_z*cos(fi)-u_y*sin(fi))/sqrt(1-u_z)+u_x*cos_theta;
                    u_y_new = sin_theta*(u_y*u_z*cos(fi)-u_x*sin(fi))/sqrt(1-u_z)+u_y*cos_theta;
                    u_z_new = -sin_theta*cos(fi)*sqrt(1-u_z)+u_z*cos_theta;
                else
                    u_x_new = sin_theta*cos(fi);
                    u_y_new = sin_theta*sin(fi);
                    if u_z >= 0
                        u_z_new = cos_theta;
                    else
                        u_z_new = -cos_theta;
                    end
                end
                u_x = u_x_new;
                u_y = u_y_new;
                u_z = u_z_new;
                % 更新电场
                sin_fi = sin(fi);
                cos_fi = cos(fi);
                % 得到A矩阵，更新m,n,s
                A = [cos_theta*cos_fi, cos_theta*sin_fi, -sin_theta;
                    -sin_fi, cos_fi, 0;
                    sin_theta*cos_fi, sin_theta*sin_fi, cos_theta];
                res = A*[m;n;s];
                [m,n,s] = deal(res(1,:), res(2,:), res(3,:));
                % 得到B矩阵，更新E1，E2
                theta_degree = ceil(rad2deg(acos(cos_theta)));
                S1 = obj.S1_vec(theta_degree);
                S2 = obj.S2_vec(theta_degree);
                FF = (abs(S2)^2*cos_fi^2 + abs(S1)^2*sin_fi^2)*(abs(E1))^2 + (abs(S2)^2*sin_fi^2+abs(S1)^2*cos_fi^2)*(abs(E2))^2 + 2*((abs(S2)^2 - abs(S1)^2))*cos_fi*sin_fi*real(E1*conj(E2));
                B = FF^(-0.5)*[S2*cos_fi, S2*sin_fi;-S1*sin_fi,S1*cos_fi];
                res = B*[E1;E2];
                [E1, E2] = deal(res(1), res(2));
                % 更新光子能量
                w = w*exp(-obj.K_a*ss);
                if scatter_type == 0
                    w = w*exp(-obj.K_sr*ss);
                else
                    w = w*exp(-obj.K_sm*ss);
                end
                % 折射次数自增
                cnt = cnt + 1;
            end
            single_context(:, cnt) = [x, y, z, E1];
            obj.context(:, :, index) = single_context;
        end
        
        % 执行函数
        function obj = run(obj)
            disp('统计光源信息...');
            R = obj.light_source.getRadius(); % 光束半径
            [X,Y] = meshgrid(linspace(-3*R, 3*R, obj.Ntr));
            [thetas, radiuses] = cart2pol(X, Y);
            total_I = 0;                      % 总光强
            for i = 1 : obj.Ntr
                disp(i);
                for j = 1 : obj.Ntr
                    x = X(j, i);
                    y = Y(j, i);
                    theta = thetas(j, i);
                    radius = radiuses(j, i);
                    I = obj.light_source.getLightIntensity(theta, radius);
                    [lambda, u_x, u_y, u_z, m, n, s, E1, E2] = obj.light_source.getInfo(theta, radius);
                    obj.light_source_info(:,j,i) = [x, y, I, lambda, u_x, u_y, u_z, m, n, s, E1, E2];
                    total_I = total_I + I;
                end
            end
            disp('开始模拟...');
            cnt = 0;
            index = 1;
            for i = 1 : obj.Ntr
                disp(i);
                for j = 1 : obj.Ntr
                    [x, y, I, lambda, u_x, u_y, u_z, mx, my, mz, nx, ny, nz, sx, sy, sz, E1, E2] = deal(obj.light_source_info(1,j,i), obj.light_source_info(2,j,i), obj.light_source_info(3,j,i), obj.light_source_info(4,j,i), obj.light_source_info(5,j,i), obj.light_source_info(6,j,i), obj.light_source_info(7,j,i), obj.light_source_info(8,j,i), obj.light_source_info(9,j,i), obj.light_source_info(10,j,i), obj.light_source_info(11,j,i), obj.light_source_info(12,j,i), obj.light_source_info(13,j,i), obj.light_source_info(14,j,i), obj.light_source_info(15,j,i), obj.light_source_info(16,j,i), obj.light_source_info(17,j,i), obj.light_source_info(18,j,i));
                    count = round(I*obj.N/total_I);
                    for k = 1 : count
                        obj.singleTimeSimulator(index, lambda, x, y, 0, u_x, u_y, u_z, [mx,my,mz], [nx,ny,nz], [sx,sy,sz], E1, E2);
                        index = index + 1;
                    end
                    cnt = cnt + count;
                end
            end
        end
        
        % 得到光强图
        function obj = showLightIntensityDiagram(obj, D, R, NNxy)
            % 统计D位置的接收点
            receiver = zeros(NNxy, NNxy);
            for i = 1 : size(obj.context, 3)
                [last_x, last_y, last_z] = deal(obj.context(1, 1, i), obj.context(2, 1, i), obj.context(3, 1, i));
                for j = 2 : size(obj.context, 2)
                    [x, y, z] = deal(obj.context(1, j, i), obj.context(2, j, i), obj.context(3, j, i));
                    if x == 0 && y == 0 && z == 0
                        break;
                    end
                    if last_z < D && z >= D   % 被无限大光屏捕获
                        % 计算与光屏交点坐标
                        x_recv = last_x + (x - last_x)*(D - last_z)/(z - last_z);
                        y_recv = last_y + (y - last_y)*(D - last_z)/(z - last_z);
                        if x_recv >= -R && x_recv <= R && y_recv >= -R && y_recv <= R % 被光屏捕获
                            % 计算交点索引
                            x_recv_index = round((x_recv + R)/(2*R/NNxy)) + 1;
                            y_recv_index = round((y_recv + R)/(2*R/NNxy)) + 1;
                            if x_recv_index > NNxy
                                x_recv_index = NNxy;
                            end
                            if y_recv_index > NNxy
                                y_recv_index = NNxy;
                            end
                            % 更新相应接收屏
                            receiver(x_recv_index, y_recv_index) = receiver(x_recv_index, y_recv_index) + 1;
                            break;
                        end
                    end
                    % 更新变量
                    [last_x, last_y, last_z] = deal(x, y, z);
                end
            end
            % 绘制光强图
            [x,y] = meshgrid(linspace(-R, R, NNxy));
            figure; 
            mesh(x*1e3, y*1e3, receiver);
            colormap(jet);
            colorbar;
            % view(2);
            % caxis([0, 300]);
            set(gca,'fontname', 'times new roman', 'fontsize',14);
            title(['{\itz}=', num2str(D), '接收处的光强'], 'fontname', '华文中宋', 'fontsize', 14);
            xlabel('\itx/mm', 'fontname', 'times new roman', 'fontsize', 14);
            ylabel('\ity/mm', 'fontname', 'times new roman', 'fontsize',14);
            zlabel('归一化强度', 'fontname', '华文中宋', 'fontsize', 14);
        end

        % 得到相位图
        function obj = showPhaseDiagram(obj, D, R, NNxy)
            k = 2*pi/obj.light_source.getLambda();
            % 统计D位置的接收点
            receiver_ph_max = zeros(NNxy, NNxy); % 用于统计接收相位
            receiver_ph_max = receiver_ph_max - pi;
            receiver_ph_min = zeros(NNxy, NNxy); % 用于统计接收相位
            receiver_ph_min = receiver_ph_min + pi;
            for i = 1 : size(obj.context, 3)
                [last_x, last_y, last_z, last_E1] = deal(obj.context(1, 1, i), obj.context(2, 1, i), obj.context(3, 1, i), obj.context(4, 1, i));
                for j = 2 : size(obj.context, 2)
                    [x, y, z, E1] = deal(obj.context(1, j, i), obj.context(2, j, i), obj.context(3, j, i), obj.context(4, j, i));
                    if x == 0 && y == 0 && z == 0
                        break;
                    end
                    if last_z < D && z >= D   % 被无限大光屏捕获
                        % 计算与光屏交点坐标
                        x_recv = last_x + (x - last_x)*(D - last_z)/(z - last_z);
                        y_recv = last_y + (y - last_y)*(D - last_z)/(z - last_z);
                        if x_recv >= -R && x_recv <= R && y_recv >= -R && y_recv <= R % 被光屏捕获
                            % 计算交点索引
                            x_recv_index = round((x_recv + R)/(2*R/NNxy)) + 1;
                            y_recv_index = round((y_recv + R)/(2*R/NNxy)) + 1;
                            if x_recv_index > NNxy
                                x_recv_index = NNxy;
                            end
                            if y_recv_index > NNxy
                                y_recv_index = NNxy;
                            end
                            % 计算交点相位
                            DD = sqrt((x_recv-last_x)^2 + (y_recv-last_y)^2 + (D-last_z)^2);
                            E1 = last_E1*exp(1i*k*DD);
                            receiver_ph_max(x_recv_index, y_recv_index) = max(receiver_ph_max(x_recv_index, y_recv_index), angle(E1));
                            receiver_ph_min(x_recv_index, y_recv_index) = min(receiver_ph_min(x_recv_index, y_recv_index), angle(E1));
                            break;
                        end
                    end
                    % 更新变量
                    [last_x, last_y, last_z, last_E1] = deal(x, y, z, E1);
                end
            end
            % 绘制最大相位图
            [x,y] = meshgrid(linspace(-R, R, NNxy));
            figure;
            mesh(x*1e3, y*1e3, receiver_ph_max);
            colormap(gray);
            colorbar;
            view(2);
            caxis([-pi,pi]);
            set(gca, 'fontname', 'times new roman', 'fontsize', 14);
            title(['{\itz}=', num2str(D), '接收处的相位（最大采样）'], 'fontname','华文中宋', 'fontsize',14);
            xlabel('\itx/mm', 'fontname', 'times new roman', 'fontsize', 14);
            ylabel('\ity/mm', 'fontname', 'times new roman', 'fontsize', 14);
            zlabel('归一化强度', 'fontname', '华文中宋', 'fontsize', 14);
            % 绘制最小相位图
            [x,y] = meshgrid(linspace(-R, R, NNxy));
            figure;
            mesh(x*1e3, y*1e3, receiver_ph_min);
            colormap(gray);
            colorbar;
            view(2);
            caxis([-pi,pi]);
            set(gca, 'fontname', 'times new roman', 'fontsize', 14);
            title(['{\itz}=', num2str(D), '接收处的相位（最小采样）'], 'fontname', '华文中宋', 'fontsize', 14);
            xlabel('\itx/mm', 'fontname', 'times new roman', 'fontsize', 14);
            ylabel('\ity/mm', 'fontname','times new roman', 'fontsize', 14);
            zlabel('归一化强度', 'fontname', '华文中宋', 'fontsize', 14);
        end

        % 统计相位串扰
        function obj = statisticalPhaseCrosstalk(obj, D, R)
            k = 2*pi/obj.light_source.getLambda();
            receiver_x_point = 1:obj.N;   % 接收点的横坐标集合
            receiver_y_point = 1:obj.N;   % 接收点的纵坐标集合
            receiver_phase = 1:obj.N;     % 接收点的相位集合
            receiver_count = 0;           % 总接收光子数
            for i = 1 : size(obj.context, 3)
                [last_x, last_y, last_z, last_E1] = deal(obj.context(1, 1, i), obj.context(2, 1, i), obj.context(3, 1, i), obj.context(4, 1, i));
                for j = 2 : size(obj.context, 2)
                    [x, y, z, E1] = deal(obj.context(1, j, i), obj.context(2, j, i), obj.context(3, j, i), obj.context(4, j, i));
                    if x == 0 && y == 0 && z == 0
                        break;
                    end
                    if last_z < D && z >= D   % 被无限大光屏捕获
                        % 计算与光屏交点坐标
                        x_recv = last_x + (x - last_x)*(D - last_z)/(z - last_z);
                        y_recv = last_y + (y - last_y)*(D - last_z)/(z - last_z);
                        if x_recv >= -R && x_recv <= R && y_recv >= -R && y_recv <= R % 被光屏捕获
                            DD = sqrt((x_recv-last_x)^2 + (y_recv-last_y)^2 + (D-last_z)^2);
                            E1 = last_E1*exp(1i*k*DD);
                            % 更新相应接收屏
                            receiver_x_point(receiver_count+1) = x_recv;
                            receiver_y_point(receiver_count+1) = y_recv;
                            receiver_phase(receiver_count+1) = angle(E1);
                            receiver_count = receiver_count + 1;
                            break;
                        end
                    end
                    % 更新变量
                    [last_x, last_y, last_z, last_E1] = deal(x, y, z, E1);
                end
            end
            % 绘制接收点相位散点图
            receiver_x_point = receiver_x_point(1:receiver_count);
            receiver_x_point = receiver_x_point .* 1000;
            receiver_y_point = receiver_y_point(1:receiver_count);
            receiver_y_point = receiver_y_point .* 1000;
            receiver_phase = receiver_phase(1:receiver_count);
            figure; 
            scatter3(receiver_x_point, receiver_y_point, receiver_phase, 1, 'k');
            % 绘制极坐标矩阵
            X = 1:receiver_count;
            Y = 1:receiver_count;
            for cnt=1:receiver_count
                [theta, ~] = cart2pol(receiver_x_point(cnt), receiver_y_point(cnt));
                theta = theta - 0.52*pi/2;
                if theta < 0
                    theta = theta + 2*pi;
                end
                X(cnt) = theta;
                Y(cnt) = receiver_phase(cnt);
            end
            Y = Y + pi;
            figure; 
            scatter(X, Y, 1, 'k');
            % 随机算法统计串扰率
            NNN = 20000000; % 随机选择2000万次
            l_1_class_cnt = 0;
            l_2_class_cnt = 0;
            l_3_class_cnt = 0;
            l_4_class_cnt = 0;
            l_5_class_cnt = 0;
            ccc = 0;
            for i = 1 : NNN
                index_0 = round(rand(1)*receiver_count);
                if index_0 == 0
                    index_0 = round(rand(1)*receiver_count);
                end
                index_1 = round(rand(1)*receiver_count);
                if index_1 == 0 || index_1 == index_0
                    index_1 = round(rand(1)*receiver_count);
                end
                ph1 = receiver_phase(index_0);
                [theta1, ~] = cart2pol(receiver_x_point(index_0), receiver_y_point(index_0));
                ph2 = receiver_phase(index_1);
                [theta2, ~] = cart2pol(receiver_x_point(index_1), receiver_y_point(index_1));
                if abs(theta1 - theta2) < 0.1 || abs(theta1 - theta2) > 0.4
                    continue;
                end
                kk = abs((ph2-ph1)/(theta2-theta1));
                % 进行归类
                if kk >= 0.5 && kk < 5.5
                    ccc = ccc + 1;
                end
                if kk >= 0.5 && kk < 1.5
                    l_1_class_cnt = l_1_class_cnt+1;
                elseif kk >= 1.5 && kk < 2.5
                    l_2_class_cnt = l_2_class_cnt+1;
                elseif kk >= 2.5 && kk < 3.5
                    l_3_class_cnt = l_3_class_cnt+1;
                elseif kk >= 3.5 && kk < 4.5
                    l_4_class_cnt = l_4_class_cnt+1;
                elseif kk >= 4.5 && kk < 5.5
                    l_5_class_cnt = l_5_class_cnt+1;
                end
            end
            disp(ccc)
            disp('l = 1: ')
            disp(l_1_class_cnt/ccc)
            disp('l = 2: ')
            disp(l_2_class_cnt/ccc)
            disp('l = 3: ')
            disp(l_3_class_cnt/ccc)
            disp('l = 4: ')
            disp(l_4_class_cnt/ccc)
            disp('l = 5: ')
            disp(l_5_class_cnt/ccc)
        end
    end
end